package com.asa.hutils.teizheng.planb;

import java.util.Arrays;

public class ASA {
	
	
	

	 public static double[] transformSecondary(double[][] householderVectors) {

	        final int m = householderVectors.length;
			double[] secondary = new double[m - 1];
			double[] main  = new double[m];

	        
	        
	        final double[] z = new double[m];
	        for (int k = 0; k < m - 1; k++) {

	            //zero-out a row and a column simultaneously
	            final double[] hK = householderVectors[k];
	            main[k] = hK[k];
	            double xNormSqr = 0;
	            for (int j = k + 1; j < m; ++j) {
	                final double c = hK[j];
	                xNormSqr += c * c;
	            }
	            final double a = (hK[k + 1] > 0) ? -Math.sqrt(xNormSqr) : Math.sqrt(xNormSqr);
	            secondary[k] = a;
	            if (a != 0.0) {
	                // apply Householder transform from left and right simultaneously

	                hK[k + 1] -= a;
	                final double beta = -1 / (a * hK[k + 1]);

	                // compute a = beta A v, where v is the Householder vector
	                // this loop is written in such a way
	                //   1) only the upper triangular part of the matrix is accessed
	                //   2) access is cache-friendly for a matrix stored in rows
	                Arrays.fill(z, k + 1, m, 0);
	                for (int i = k + 1; i < m; ++i) {
	                    final double[] hI = householderVectors[i];
	                    final double hKI = hK[i];
	                    double zI = hI[i] * hKI;
	                    for (int j = i + 1; j < m; ++j) {
	                        final double hIJ = hI[j];
	                        zI   += hIJ * hK[j];
	                        z[j] += hIJ * hKI;
	                    }
	                    z[i] = beta * (z[i] + zI);
	                }

	                // compute gamma = beta vT z / 2
	                double gamma = 0;
	                for (int i = k + 1; i < m; ++i) {
	                    gamma += z[i] * hK[i];
	                }
	                gamma *= beta / 2;

	                // compute z = z - gamma v
	                for (int i = k + 1; i < m; ++i) {
	                    z[i] -= gamma * hK[i];
	                }

	                // update matrix: A = A - v zT - z vT
	                // only the upper triangular part of the matrix is updated
	                for (int i = k + 1; i < m; ++i) {
	                    final double[] hI = householderVectors[i];
	                    for (int j = i; j < m; ++j) {
	                        hI[j] -= hK[i] * z[j] + z[i] * hK[j];
	                    }
	                }
	            }
	        }
	        main[m - 1] = householderVectors[m - 1][m - 1];
	        
	        
	        return secondary;
	    }
	
	
	
	
	
	/**
	 * Q矩阵
	 * @param householderVectors
	 * @return
	 */
	public double[][] getQT(double[][] householderVectors) {
           final int m = householderVectors.length;
           double[][] qta = new double[m][m];
           double[] secondary = transformSecondary(householderVectors);
           // build up first part of the matrix by applying Householder transforms
           for (int k = m - 1; k >= 1; --k) {
               final double[] hK = householderVectors[k - 1];
               qta[k][k] = 1;
               if (hK[k] != 0.0) {
                   final double inv = 1.0 / (secondary[k - 1] * hK[k]);
                   double beta = 1.0 / secondary[k - 1];
                   qta[k][k] = 1 + beta * hK[k];
                   for (int i = k + 1; i < m; ++i) {
                       qta[k][i] = beta * hK[i];
                   }
                   for (int j = k + 1; j < m; ++j) {
                       beta = 0;
                       for (int i = k + 1; i < m; ++i) {
                           beta += qta[j][i] * hK[i];
                       }
                       beta *= inv;
                       qta[j][k] = beta * hK[k];
                       for (int i = k + 1; i < m; ++i) {
                           qta[j][i] += beta * hK[i];
                       }
                   }
               }
           }
           qta[0][0] = 1;
   
       // return the cached matrix
       return qta;
   }
	
	
	

}
